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Chiral p-wave superconducting state supports a rich spectrum of topological excitations different 
from those in conventional superconducting states. Besides domain walls separating different chiral 
states, chiral p- wave state supports both singular and coreless vortices also interpreted as skyrmions. 
Here, we present a numerical study of the energetic properties of isolated singular and coreless 
vortex states as functions of anisotropy and magnetic field penetration length. In a given chiral 
state, single quantum vortices with opposite winding have different energies and thus only one 
kind is energetically favoured. We find that with the appropriate sign of the phase winding, two- 
quanta (coreless) vortices are always energetically preferred over two isolated single quanta (singular) 
vortices. We also report solutions carrying more flux quanta. However those are typically more 
energetically expensive/metastable as compared to those carrying two flux quanta. 


Chiral p -wave superconducting state is an exotic state 
that, in addition to usual U(l) gauge symmetry, spon¬ 
taneously breaks time-reversal symmetry. Higher bro¬ 
ken symmetries there, implies a much richer spectrum of 
topological excitations as compared to conventional su¬ 
perconducting and superfluid states. Chiral p -wave pair¬ 
ing is realized in the A-phase of superfluid 3 He, were 
variety of complex topological defects were investigated 
[1-8]. In a superconducting p-wave state, due to the cou¬ 
pling to the vector potential, topological defects exhibit 
different properties. This coupling affects their energy 
and determines their role in the magnetic properties of 
such superconductors. Layered perovskite superconduc¬ 
tor Sr2Ru04 is a candidate material where various exper¬ 
imental evidences suggest possible realization of a p -wave 
superconducting state [9, 10]. Similar models were also 
considered in connection to the superconducting state of 
heavy fermion compound UPt3 [11, 12] (see e.g. [13, 14] 
for recent discussion of superconducting state in that ma¬ 
terial). 

Spontaneous breaking of time-reversal symmetry for 
chiral p -wave state, implies the existence domain walls 
that separate regions with two different time-reversal 
symmetry broken (TRSB) ground-states. These domain 
walls, support spontaneous supercurrent that can gener¬ 
ate magnetic fields [15-20]. The domain walls in chiral p- 
wave superconductors could be created via Kibble-Zurek 
mechanism, and these properties can be used for their 
control [21]. However, in Si'2Ru04 no indication of such 
a held was found in magnetic imaging microscopy ex¬ 
periments [22-24]. Besides domain walls, chiral p-wave 
superconductors feature rich spectrum of topological de¬ 
fects including various vortices and skyrmions [25, 26]. 
Related topological defects were also discussed in the con¬ 
text of heavy fermion superconductor UPt3 [27]. 

In zero held, both chiral (ground-)states are degener¬ 
ate in energy and this degeneracy is lifted by an exter¬ 
nally applied magnetic held along the c-axis. For a given 
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sign of the magnetic held parallel to the c-axis, only one 
of the chiral states is stable and the time-reversed state 
is energetically penalized. Likewise, vorticity of the su¬ 
perconducting condensates lifts the degeneracy between 
both chiral states. When the dominant component forms 
a vortex, it induces the time-reversed (subdominant) chi¬ 
ral component, in the vicinity of the core. The winding 
of the induced component is not independent of that of 
the dominant component. It has a 47 t winding of the rel¬ 
ative phase between components, that follows from the 
Cooper pairs having nonzero internal orbital momentum 
[28]. Since the magnetic held lifts the degeneracy be¬ 
tween chiralities, vortices and anti-vortices have different 
properties [29]. 

It is experimentally seen that in an applied external 
held, Sr2Ru04 exhibits vortex lattices with square sym¬ 
metry at high helds [30-32], and a transition to trian¬ 
gular lattice in lower helds [32, 33]. Earlier theoretical 
calculations based on Ginzburg-Landau model for chiral 
p -wave superconductivity in Sr2Ru04 [34-36], are consis¬ 
tent with these observed transitions of the vortex lattice 
structure. 

Besides single-quanta vortices, there also exists vor¬ 
tices carrying multiple quanta of the magnetic flux and 
that, as they are coreless, are essentially different from 
single-quanta vortices. For example as discussed in more 
details below, the component induced by a doubly quan¬ 
tized vortex in the dominant component has zero wind¬ 
ing in subdominant one [29]. In this paper we demon¬ 
strate that the two-quanta (coreless) vortices, which can 
also be denoted as skyrmions, are energetically favoured 
as compared to (isolated) single-quanta vortices. Ear¬ 
lier works in the context of UPt 3 , even claim that lat¬ 
tices of similar two-quanta vortices may be energetically 
favoured as compared to those of single quanta [37, 38]. 
The possible existence of lattices of different coreless vor¬ 
tices carrying single flux quantum in UPt3 was also dis¬ 
cussed recently [14]. It was also recently shown in the 
context of S^RuCL, based on solutions of microscopic 
Eilenberger equations, that lattices of two-quanta vor¬ 
tices are favoured for certain parameter sets [39]. Yet, 
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such lattices of two-quanta vortices were never observed 
in Sr 2 Ru 04 . This motivates this work to further investi¬ 
gate the thermodynamic stability of skyrmions for broad 
parameter range. 

In a previous work [40] , we reported isolated skyrmion 
solutions in a model for chiral p -wave superconductor. 
For the studied case of one of the chiralities, skyrmions 
can be energetically favoured as compared to vortices 
(see also remark [41]). The skyrmions carrying two flux 
quanta are directly related to the two-quanta vortices 
mentioned above. However it was also demonstrated in 
Ref. [40] that there are (meta-)stable skyrmions carry¬ 
ing larger number of flux quanta. In this model the en¬ 
ergy and structure of vortices and skyrmions depends 
on the chirality. Equivalently, for a given chiral state, 
vortex/skyrmion solutions are not the same as anti- 
vortex/anti-skyrmion. It thus calls for further investi¬ 
gation of vortex and skyrmion solutions (carrying two 
and more quanta), which we present below. 

In the coordinate system where the crystal anisotropy 
axis is c || z, the order parameter of the p x + ip y 
state is described by a complex two-dimensional vector 
V = {Vx,Vy) [10, 11, 42]. Introducing the chiral order 
parameter components = (r] x ± irj v ) /a/ 2 , the dimen¬ 
sionless Ginzburg-Landau free energy reads as (see e.g. 
[34-36]): 

F = \V x A\ 2 + \D v+ \ 2 + \D V _\ 2 (la) 

+ iy + l)Re [(D x r)+)*D x r]_ - (D y r)+)*D y r)-] (lb) 

+ {y- l)Im l(D x i 1+ yD v ri- + (D y r) + )*D x r)-\ (lc) 

+ 2\rj + rj-\ 2 + pR e(r]fr] 2 _) + ^ -|? 7 a | 2 + *|? 7 a| 4 - (Id) 

a=± 

There we use dimensionless units were the free energy is 
normalized to the condensation energy, and the lengths 

_I /o 

are given in units of £ = (op (T — T c )) ' . The mag¬ 
netic field B = V x A is in units of B c = <bo/(27rA£). 
In these units, the gauge coupling g that appears in 
the covariant derivative D = V + igA is related to 
the Ginzburg-Landau parameter g~ x = A/£. The free 
energy (1) was derived from the weak coupling micro¬ 
scopic theory [34, 35]. The anisotropy parameter v de¬ 
termines the anisotropy in the iry-plane (|^| < 1 for 
the energy to be positively defined). It measures the 
tetragonal distortions of the Fermi surface, which has 
cylindrical geometry for v = 0, and is defined as v = 
i(v x ) - 3(v 2 v 2 ))/((vl) + (v 2 x v 2 y )) (where (•) denote aver¬ 
age over the Fermi surface). In the model Eq. (1), the 
dependence on the third coordinate is not considered (i.e. 
assuming two-dimensional system or translational invari¬ 
ance along 0 -axis). 

The ground-state that minimizes the potential terms 
in ( 1 ) is degenerate and the solutions are (? 7 +,? 7 _) = 
(1,0) and (0,1). The theory (1) is invariant under the 

(discrete) time-reversal operations T, as {rj±,B} —y 
{ 77 ^-, — B}. This invariance is spontaneously broken by 
the ground-state. The spontaneous breakdown of the 


discrete time-reversal symmetry dictates that the theory 
allows domain wall solutions that interpolate between re¬ 
gions with different ground-states. Such domain walls, 
rather generically created during phase transition where 
the discrete symmetry is broken [ 21 ], carry a magnetic 
field perpendicular to the ccy-plane [17, 18]. The discrete 
degeneracy of the ground state is lifted by the magnetic 
field. Thus, depending on the direction of the external 
field, only one state is stable. Likewise, the vorticity of 
the superconducting condensates lifts the degeneracy be¬ 
tween chiral (ground-)states. 

As the components 77 + and rj- behave differently for 
different sign of the winding, a complete study requires 
to consider both situations of counter-clockwise (pos¬ 
itive) and clockwise (negative) vorticities. Note that 
this is equivalent to considering only positive vorticity 
but for both chiral states. For example, the configu¬ 
ration with a winding n+ = +1 on the ground-state 
(? 7 +,? 7 _) = ( 1 , 0 ) can be obtained by applying the time- 
reversal operation T on the configuration whose ground 
state is (rj+,r]-) = (0,1) with the winding n_ = —1. In 
the following, we choose to fix the dominant component 
of the order parameter to be rj- [i.e. the ground state 
is (? 7 + ,? 7 _) = ( 0 , 1 )] and thus need to investigate both 
positive and negative vorticities. 

The asymptotic vorticity of the dominant component 
? 7 _ determines the sign of B z , as well as the vorticity of 
the subdominant component 77 + [35], according to: 

? 7 _ oc e m ~ 0 , 77 + oc e m+e and n + =n-+2, ( 2 ) 

where 9 is the polar angle. The relative phase between g + 
and Tj— ( 2 ), follows from the internal orbital momentum 
of Cooper pairs. In the Ginzburg-Landau model (1), this 
is the structure of mixed gradient that constraints the 
relative phase. Note that since the subdominant com¬ 
ponent r)+ vanishes asymptotically [i.e. it recovers its 
ground state value rj+ = 0 ], the winding n+ can be lo¬ 
cated only in a close vicinity of the vortex core. Note also 
that the winding of the subdominant component does not 
affect the overall flux quantization, because the density 
of that component vanishes away from the vortex. From 
( 2 ) it is rather straightforward to see that the vortex with 
the vorticity (n+,n-) = (+3,+l) and the (anti-)vortex 
with (tt + ,71_) = (+1,-1) have different core structures 
and thus different energy. 

In order to investigate the energetic properties of vor¬ 
tex matter, the fields r]± and A are discretized using 
a finite-element framework [43] and the free energy (1) 
is minimized using an nonlinear conjugate gradient al¬ 
gorithm (see the appendix for details). In simulations 
of chiral p-wave superconductors on a finite domain, a 
special attention is required for boundary conditions in 
order to yield edge currents (see for example discussions 
in [19, 42, 44]). Here, we are interested in the intrinsic 
energetic properties of isolated defects. Thus vortex con¬ 
figuration is created by an initial guess and placed them 
at the center of a large domain, with open boundary con¬ 
ditions, letting the fields freely reach the ground-state. 
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Figure 1. (Color online) - Vortex states for the windings (n+,n_) of the components 77 + and Note that the winding of the 
dominant component (here g~), specifies the flux carried by the vortex configuration. The parameters of the Ginzburg-Landau 
functional are g = 0.3 and v — 0.2. The first line shows the magnetic field B, while second and third line respectively display 
|? 7 _| and |? 7 +|. The fourth line shows the relative phase i p- — between and g-. Winding of the relative phase indicates 
the position of the cores of g + and ??_. The first block shows vortex solutions carrying one to four flux quanta with B z > 0, 
while the second block shows the corresponding vortices with B z < 0. 


By choosing a sufficiently large domain, this ensures that 
within numerical accuracy vortices will not interact with 
boundaries and thus we are able to probe their intrin¬ 
sic structure and energy properties, without effects of 
boundary conditions. As it specifies the topological sec¬ 
tor, a starting configuration with a given winding ?z_ 
of the dominant component rj- leads, after convergence 
of the algorithm, to a configuration that behaves as ex¬ 
pected from (2). We systematically construct vortex so¬ 
lutions carrying one to four flux quanta for parameter 
space defined by wide range of values of the g and u. 
Fig. 1 shows typical vortex solutions with different vor- 
ticities. Along this paper we also refer to vortices carry¬ 
ing multiple flux quanta, as skyrmions. The reason for 
that terminology is that they have additional topological 
properties, as compared to single quanta vortices. This 
is explained in more details by the end of the paper. 

The first and second blocks in Fig. 1 respectively show 
vortex solutions with B z > 0 and B z < 0. Vortices car¬ 
rying from one to four flux quanta are displayed within 
each block. As expected from Eq. (2), single winding of 
the dominant component induces core structure of the 
subdominant component with different winding depend¬ 
ing on that of ??_ (see the first column of each block). It is 
instructive to consider the last row in Fig. 1, that displays 
the relative phase between and g + . In agreement with 
( 2 ), asymptotically the relative phase p- — ip + = —2 9, 
reflecting the orbital angular momentum difference be¬ 
tween rj _ and g+. Moreover, the relative phase also in¬ 
dicates the position of the singularities of the compo¬ 
nents of the order parameter. Remarkably single quanta 


vortices are singular defects, since singularities in both 
components overlap (and thus g+ = g~ = 0). On the 
other hand, since both components never simultaneously 
vanish, two-quanta vortices are coreless defects. Inter¬ 
estingly the n_ = —2 configuration features a ir jump of 
the relative phase when going outward from the vortex 
core. Inside the vortex core the time-reversed chiral state 
(ry+,? 7 _) = ( 1 , 0 ) is induced, while the (g + ,g_) = ( 0 , 1 ) 
state is recovered asymptotically. The two quanta vor¬ 
tices thus feature a domain wall when going away from 
the center. The 7 r jump of the relative phase for the 
n_ = —2 is located at this domain wall. 

Like in conventional superconductors, the magnetic 
field for single quanta vortices is localized at the vor¬ 
tex core and screened at length scales determined by the 
penetration depth A. Interestingly, the magnetic field for 
two-quanta vortices, and especially for n_ = — 2 , is not 
homogeneously distributed in the core. Rather it is local¬ 
ized at a given distance from the center and spread along 
the ring of the domain wall. Note that similar vortex 
configurations were also found to exist in the context of 
two-component model with U(l) x U(l) x Z 2 symmetry 
[45]. The ring-like distribution of the magnetic field for 
the two-quantum vortex can be understood as follows: 
B outside the vortex is screened by the (partial) cur¬ 
rents in g_ that run counter-clockwise, while inside the 
vortex currents in g + are responsible for the screening. 
Since g+ vanishes away from the vortex core, it cannot 
contribute to the screening asymptotically. Conversely, 
? 7 _ vanishes at the vortex core and this is the induced 
subdominant component 77 + that screens B close to the 
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Figure 2. (Color online) - Left panel shows the ratio of the energies E{n- = —1) of n_ = —1 single quanta vortices and 
E(n- = +1) of n_ = +1 vortices, as functions of the anisotropy parameter v and of the gauge coupling g. This ratio is 
always smaller than 1 implying that (n+,n_) = (+1,-1) are always less energetic than (n+,n_) = (+3,+1). The right panel 
displays the ratio of the energies 2 E(n- = —1) of two n_ = —1 single-quanta vortices and E(n- = —2) of one n_ = —2 
vortex, as functions of the anisotropy parameter v and of the gauge coupling g. This ratio is always larger than 1 implying 
that two-quanta vortices are less energetic than two isolated single-quanta vortices. 


center of the vortex core. The reason it can contribute to 
screening (inside the vortex) without having vorticity on 
its own, is only due to supercurrent produced by the vec¬ 
tor potential (like the Meissner currents on the bound¬ 
ary of ordinary superconductors). Since those currents 
circulate clockwise, they compensate with the currents 
in r]_ so that at a certain distance (at the domain wall) 
there is no screening current. The magnetic field is thus 
localized at the domain wall. Although the core struc¬ 
ture of single-quanta vortices are different depending on 
the sign of n_, their profile of the magnetic field looks 
quite similar. When considering vortices with |n_| > 1, 
both the core structure and the magnetic field profile are 
strikingly different and the skyrmions with negative n_ 
do not resemble those with positive n_. Apart from the 
n_ = —2, —3 skyrmions, the configurations that carry 
multiple flux quanta are far from being axially symmet¬ 
ric. Note that the n_ = —4 skyrmions resembles as some 
kind of bound state of two n _ = —2 skyrmions. As 
we will discuss below, this makes their decay into two 
n_ = —2 vortices rather easy. 

Since the core structure is different, it is quite natural 
to expect that, unlike in conventional superconductors, 
vortices with opposite winding (and thus opposite direc¬ 
tions of the magnetic field) are non-degenerate in energy. 
The (n + ,n_) = (+3,+1) vortex has more total vorticity 
than the (n+,n_) = (+1,-1). Thus one could naively 
expect that the n_ = — 1 vortices would be favoured as 
compared to ?r_ = +1. We systematically compared the 
energies of both single-quanta vortices for all values of 
the anisotropy parameter v and of the gauge coupling 
g. The diagram in Fig. 2(a), shows that the ratio of the 
energies of the single quanta vortices with n_ = — 1 and 
n_ = +1, is always less than one. This implies that 
vortices = — 1 are always energetically favoured, as 
compared to those with n_ = +1. The first critical field 
of a vortex carrying a flux 4> is H c i = E/ 2$, where E 
is its energy. As a result, Fig. 2(a) also implies that the 
n_ = — 1 vortices also have lower first critical field H c i in 


agreement with Refs. [36, 46]. Although both n_ = ±1 
are perturbatively stable (i.e. they are minima of J -), 
only n_ = — 1 is absolutely stable. 

Note that the naive estimates based on counting the to¬ 
tal vorticity provide the correct picture that (n+,n_) = 
(+1,-1) vortices are less energetic than (n+,n_) = 
(+3, +1) ones. It thus makes sense to apply the same 
arguments to configurations carrying more than one flux 
quantum. In the sector with negative n_, there are 
two possibilities to make a configuration that carries two 
flux quanta. Either to create two isolated (n + , n_) = 
(+1,-1) vortices carrying one flux quantum each or 
to create one (n+,n_) = (0,-2) two-quantum vortex. 
It turns out that a two-quantum vortex with smaller 
number of singularities is favoured as compared to two 
isolated single-quanta vortices. Fig. 2(b) displays the 
ratio of the energies of two (isolated) n— = — 1 vor¬ 
tices and one n_ = —2 vortex. This ratio is always 
larger than one, thus implying that two-quanta vortices 
are energetically favoured as compared to two isolated 
single-quanta vortices. Note that the quantity displayed 
in Fig. 2(b), is actually also the ratio of first critical 
fields associated with single and double quanta vortices 
H c i(n_ = —1)/H c i(n_ = —2). Note also that smaller 
H c i for a higher-flux vortex does not necessarily imply 
that such vortices will nucleate first in low magnetic field. 
That is, due to higher winding they carry larger mag¬ 
netic flux and thus can have a higher potential barrier to 
enter the sample (compared with the discussion of Bean- 
Livingston barrier in single component superconductors 
[47]). The vortices (n_ = —1) and (n_ = —2) should 
interact differently with the Meissner currents and im¬ 
age charges, and thus even if the (n_ = —2) vortices 
have lower H cl , the interaction with the boundary may 
instead favour the entry of the vortices with (n_ = — 1). 

We also calculated the energy diagram similar to that 
in Fig. 2(b), but for vortices carrying three flux quanta 
n_ = —3 (data not shown). We found that unlike for 
n_ = —2, the n_ = —3 are not always stable. That 
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Figure 3. (Color online) - Ratio of the energies of multiple quanta vortices with n_ > 0 compared to isolated vortices carrying 
smaller number of flux quanta, as functions of the anisotropy parameter v and of the gauge coupling g. The solid line separates 
the regions where isolated vortices with smaller flux are preferred over single vortex carrying more flux quanta. Panel (a) and 
(b) respectively show the energetic behaviour of double (resp. triple) quantum vortex compared to isolated vortices with single 
flux quantum. Below the solid line, that is for more isotropic (small \u\) or more type-2 (small g), isolated single quanta vortices 
are energetically favoured. Panel (c) shows the energy ratio of two (isolated) triple-quanta vortices compared to three double 
quanta. This indicates subtle sublinear energy scaling with the number of flux quanta. 


is, in some regions of the parameter space the = —3 
is found, but in some other regions it decays into one 
single-quantum plus one double-quantum vortex. We 
find that for n_ < 0, the = —2 skyrmions are all 
energetically favoured. This behaviour can already be 
anticipated from the last column in Fig. 1 where the four 
quanta n_ = —4 skyrmion seems to be a bound state of 
two ri- = —2 vortices. As the skyrmions with < —2 
are more energetic than those with n_ = —2, one can 
easily see that the n_ = —4 configuration can decay into 
two n_ = —2 vortices and thus reduce its total energy. 
We find that in the regions where the n_ = — 3 vortices 
exist, they are more energetic than three isolated sin¬ 
gle quanta vortices (or one single plus one two-quantum 
vortex). 

Although being energetically unfavoured, it is still in¬ 
structive to consider the properties of vortices carrying 
multiple flux quanta, with n _ > 0. Diagrams in Fig. 3 
show the ratio of the energies of multiple quanta vor¬ 
tices with n_ > 0, compared to that of isolated vortices 
carrying smaller flux. The situation for > 0 is ac¬ 
tually very different from that with n_ < 0. Panels (a) 
and (b) in Fig. 3 display the ratio of the energies of iso¬ 
lated single-quanta vortices with that of vortices carrying 
two and three flux quanta. Depending on the anisotropy 
parameter v and on the gauge coupling g this ratio can 
either be smaller or larger than one and the solid lines on 
the diagram show when these are degenerate in energy. 
Below the solid line, isolated single-quanta vortices are 
energetically favoured as compared to multi-quanta vor¬ 
tices. Above this line, these are the vortices carrying two 
or three flux quanta which are favoured. Thus tetrago¬ 
nal distortions of the Fermi surface (i.e. larger \v\) tend 
to favour n_ = +2 (and to a lesser extend n_ = +3), 
as compared to isolated ?r_ = +1 vortices. Note that 
the solid lines in panel (a) and (b) do not coincide. The 
panel (c) shows the comparison between three isolated 
double quanta and two isolated triple quanta vortices. 
Here again, depending on v and g 1 either can be pre¬ 
ferred. This suggests complicated sublinear scaling of 


the energy with the number of flux quanta. 

The coreless nature of the two-quanta vortices implies 
that these have additional topological properties that are 
absent for single-quanta vortices. If the order parame¬ 
ter fj = (? 7 + ,? 7 _) does not vanish (fj / 0), a pseudo¬ 
spin (unit) vector n can be defined as the projection 
of the order parameter on spin-1/2 Pauli matrices cr: 
m = tfaifj/rffj (see detailed discussion of the pseudo¬ 
spin formalism for multi-component GL models in [48]). 
Fig. 4 shows pseudo-spin texture for vortex solutions cor¬ 
responding to those displayed in Fig. 1. Note that n is ill- 
defined for singular vortices, since there fj = 0 at the core 
(i.e. singularities in both components overlap). Coreless 
vortices on the other hand have well defined pseudo-spin 
projection which is a map n : R 2 —► S 2 . Since at spa¬ 
tial infinity, n = (0,0,—1), the plane R 2 can be com- 
pactified to S 2 so that the pseudo-spin becomes a map 
n : S 2 — >• S 2 . The homotopy invariants ^(S 2 ) £ TL as¬ 
sociated with such maps defines the integer-valued topo¬ 
logical charge 

Q = — [ n ■ d x n x d y n dx dy , (3) 

w J R2 

which can be used to classify various field configurations. 
Heuristically, Q counts the number of times that the tar¬ 
get sphere S 2 is wrapped while covering the xy-plane. 
Singular configurations for which the pseudo-spin is not 
everywhere well-defined, have Q = 0. Non-singular solu¬ 
tions on the other hand, and in particular coreless vor¬ 
tices, have Q = g<&/2-K £ 7L (where $ is the flux). For 
example the two-quanta vortices, which are coreless, are 
characterized by Q = 2. The fact that Q = 0 for singular 
vortices can easily be seen from the plot of the pseudo¬ 
spin texture Fig. 4. There n never reaches the north pole 
and thus do not fully cover the unit sphere. 

Here, we reported a large-scale numerical investigation 
of the energy properties of isolated single and multiple 
quanta vortices/skyrmions in a Ginzburg-Landau model 
of chiral p- wave superconducting state. As pointed out 
previously, for a given ground-state chirality, vortices and 
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Figure 4. (Color online) - Pseudo-spin texture n defined as the projection of superconducting degrees of freedom onto spin-1/2 
Pauli matrices. Different panels correspond to the solutions with different vorticity displayed in Fig. 1. First line shows vortices 
with B z > 0 and the second line is for B z < 0. The multi-quanta skyrmions are characterized by the topological charge (3) 
Q = n_. Single quanta vortices on the other hand, do not cover the whole sphere (i.e n z < n“ ax < 1) thus they have vanishing 
charge Q = 0. 


anti-vortices are inequivalent. Thus we performed study 
for both orientations of the winding. The vortices with 
winding ro_ = — 1 in the dominant component are always 
preferred to those with winding n_ = +1. We also found 
that vortices carrying two flux quanta with n_ = — 2 are 
always energetically favoured as compared to two iso¬ 
lated single-quanta vortices. Vortices with higher flux 
and negative n_, on the other hand, are either unsta¬ 
ble or have higher energies per flux quantum. We also 
reported the structural and energetic properties of (meta- 
)stable skyrmions with various topological charge (i.e. for 
n_ > +1). The calculations show complicated sublinear 
scaling of the energy with the number of flux quanta that 
qualitatively agrees with previous works for a smaller pa¬ 
rameter set in a related model [40]. Due to their very 
characteristic profiles of the magnetic field, their exper¬ 
imental observation, in e.g. scanning Hall and scanning 
SQUID experiments would provide a strong evidence of 
chiral p -wave superconductivity in the candidate materi¬ 
als described by the model (1). Note however that various 
aspects of microscopic physics may alter the form of the 
Ginzburg-Landau model (1), and in particular the bal¬ 
ance between the different coefficients entering the free 
energy. This is currently a subject of ongoing studies, 
in connection with Sr 2 Ru 04 (see e.g. [19, 20]). Note 
added: after the completion of this work a study ap¬ 
peared reporting stable skyrmions as well as vortices, in 
this model, affected by mesoscopic effects in small sam¬ 
ples [49]. 
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Appendix A: Numerical Methods 

In this work we used the dimensionless two-component 
Ginzburg-Landau theory Eq.(l) that was previously mi¬ 
croscopically derived in the weak coupling limit (see for 
example Refs. [34, 35]). In this work, we focus on the 
properties of vortex solutions in the a;y-plane and neglect 
the dependence over the third coordinate z. This means 
that our solutions are either purely two dimensional, or 
describe bulk configuration, assuming translation invari¬ 
ance along x-axis (and thus neglecting possible surface 
effects). 

For the numerical investigation, the two-dimensional 
problem (1) is defined on a bounded domain Q C R 2 . 
The boundary conditions for chiral p-wave superconduc¬ 
tors can be very involved. Namely, in order to simulate 
chiral p -wave superconductors on a finite domain, a spe¬ 
cial attention has to be paid to boundary conditions to 
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take into account edge currents properties. However, we 
are interested here in the intrinsic energetic properties 
of isolated defects. Thus we consider isolated vortices 
in large grids (such that there are no interactions with 
boundaries) and let the fields freely recover the ground- 
state. As a result, we probe the intrinsic structure and 
energy properties of vortices without any deformation 
originating from boundary behaviour. The simulation is 
run for a zero applied field (so that there are no Meissner 
currents), and the flux carrying solution is generated by a 
starting condition with a given winding of the dominant 
component. Because it enjoys topological protection, the 
(dominant) component cannot unwind by means of con¬ 
tinuous transformations and thus topological properties 
(winding of the dominant component) are preserved by 
an energy minimization algorithm. Note that as simula¬ 
tions are run on a large but finite domain, there is still 
a possibility to change the topological sector, by moving 
the vortex across the boundary. This is possible, because 
without external fields there are no Meissner currents to 
prevent escape of a vortex. Note however that as we 
choose to work with large grids, the vortices in practice 
do not interact with boundaries, and thus they do not 
escape from the domain. The advantage of this choice is 
that it is guaranteed that obtained solutions are not af¬ 
fected by boundaries and that the calculated energies are 
those of isolated defects. The configurations displayed in 
the paper are close-up views of these defects. 

For the actual numerical computation, the variational 
problem of minimizing the free energy is defined using 
a finite element formulation provided by the Freefem++ 
library [43]. Discretization within finite element formu¬ 
lation is done via a (homogeneous) triangulation over 
Q. based on Delaunay-Voronoi algorithm. Functions are 
decomposed over a continuous piecewise quadratic ba¬ 
sis on each triangle. The accuracy of such method is 
controlled through the number of triangles, (we typically 
used 3 ~ 6 x 10 1 2 * 4 5 6 7 ), the order of expansion of the basis 
on each triangle ( 2 nd order polynomial basis on each tri¬ 
angle), and also the order of the quadrature formula for 
the integral on the triangles. A nonlinear conjugate gra¬ 
dient algorithm is used to solve the variational nonlinear 


problem (i.e. to find the minima of T). The algorithm 
is iterated until relative variation of the norm of the gra¬ 
dient of the functional T with respect to all degrees of 
freedom is less than 10~ 8 (we verified that for this value, 
the configuration does not evolve and the energy remains 
constant). 

For the minimization procedure to lead to a configu¬ 
ration that has the expected topological properties, the 
starting field configuration should exhibit itself those de¬ 
sired topological properties. Although strictly speak¬ 
ing there is no infinite energy barrier between different 
topological sectors in finite domains, the barrier is finite 
but large enough to prevent any unwinding. Thus typ¬ 
ically gradient minimization converges to the configura¬ 
tion that has the topological properties of the starting 
guess. In order to have efficient numerics, it is also im¬ 
portant that the starting field configuration is the closest 
as possible to the minimal energy configuration. The 
initial field configuration carrying N v flux quanta is pre¬ 
pared by using an ansatz which imposes phase winding 
of the dominant component (??_ = | 77 _|e lv_ ) around a 
given point (x k ,y k ): 


\V -1 = l| (l +tanh (J-{n k (x,y) - (Al) 

N v 

<p_ = tan -1 

k -1 

where Tl k (x,y) = y/(x- x k ) 2 + (y - y k ) 2 and £ a 

parametrizes the core size. The parametrization of 77 +, 
with nonzero density in the core enhances the conver¬ 
gence to form coreless defects. Finally, the starting 
configuration for the vector potential of the magnetic 
field A , is determined by solving Ampere’s law equation 
V x B + J = 0, for the supercurrent J = ST/SA spec¬ 
ified by the superconducting condensates given by (Al). 
Being an equation linear in A , this operation is rapidly 
solved. Once the starting configuration is constructed, 
all degrees of freedom are relaxed simultaneously. 


y~Vk 
x - x k 


and | ? 7 + 1 = 1 - |t?_| , (A2) 
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